---
title: 'Between ideology and security interests: the influence of political parties
  on military missions in Germany'
author: "Christina Stremming"
output: html_document
---

```{r}
# Load packages
library(dplyr)
library(readr)
library(car)
library(ResourceSelection)
library(sandwich)
library(lmtest)
library(pscl)
library(pROC)
library(knitr)
library(ggplot2)
library(broom)
library(coefplot)
library(xtable)
library(stargazer)
library(knitr)
library(kableExtra)


setwd("C:/Users/Stremming/Documents/Paper 1/")
load("firstdraft.Rda")
``` 

```{r data-preparation}
#Factor variable for IO missions
missionvote$NATO <- relevel(as.factor(missionvote$NATO), ref = "8")

# Table 2 Appendix Distribution of missions
print(table(missionvote$support_dummy))
print(table(missionvote$NATO))

# Table 4 Appendix Voting behavior across parties
count_votes <- function(mission_type, party_data) {
  yes_votes <- sum(party_data$support_dummy == 1 & party_data[[mission_type]] == 1)
  no_votes <- sum(party_data$support_dummy == 0 & party_data[[mission_type]] == 1)
  return(paste(yes_votes, "/", no_votes))
}

table_votes <- data.frame(
  Organisation = c("NATO", "UN", "EU", "Others"),
  AfD = c(count_votes("NATO_dummy", subset(missionvote, party == "0")),
          count_votes("UN_dummy", subset(missionvote, party == "0")),
          count_votes("EU_dummy", subset(missionvote, party == "0")),
          count_votes("Other_dummy", subset(missionvote, party == "0"))),
  Greens = c(count_votes("NATO_dummy", subset(missionvote, party == "3")),
             count_votes("UN_dummy", subset(missionvote, party == "3")),
             count_votes("EU_dummy", subset(missionvote, party == "3")),
             count_votes("Other_dummy", subset(missionvote, party == "3"))),
  CDU_CSU = c(count_votes("NATO_dummy", subset(missionvote, party == "4")),
              count_votes("UN_dummy", subset(missionvote, party == "4")),
              count_votes("EU_dummy", subset(missionvote, party == "4")),
              count_votes("Other_dummy", subset(missionvote, party == "4"))),
  Left = c(count_votes("NATO_dummy", subset(missionvote, party == "6")),
           count_votes("UN_dummy", subset(missionvote, party == "6")),
           count_votes("EU_dummy", subset(missionvote, party == "6")),
           count_votes("Other_dummy", subset(missionvote, party == "6"))),
  FDP = c(count_votes("NATO_dummy", subset(missionvote, party == "13")),
          count_votes("UN_dummy", subset(missionvote, party == "13")),
          count_votes("EU_dummy", subset(missionvote, party == "13")),
          count_votes("Other_dummy", subset(missionvote, party == "13"))),
  SPD = c(count_votes("NATO_dummy", subset(missionvote, party == "23")),
          count_votes("UN_dummy", subset(missionvote, party == "23")),
          count_votes("EU_dummy", subset(missionvote, party == "23")),
          count_votes("Other_dummy", subset(missionvote, party == "23")))
)

kable(table_votes, format = "html", caption = "Cumulative Voting Results by Parties and International Organizations") %>%
  kable_styling(full_width = F)


#Tabelle 5 Appendix Abstimmungverhalten Parteien intern 
missionvote$yes.votes <- as.numeric(missionvote$yes.votes)
missionvote$no.votes <- as.numeric(missionvote$no.votes)

count_votes <- function(mission_type, party_data) {
  yes_votes <- sum(party_data$yes.votes[party_data[[mission_type]] == 1])
  no_votes <- sum(party_data$no.votes[party_data[[mission_type]] == 1])
  return(paste(yes_votes, "/", no_votes))
}

table_votes <- data.frame(
  Organisation = c("NATO", "UN", "EU", "Others"),
  AfD = c(count_votes("NATO_dummy", subset(missionvote, party == "0")),
          count_votes("UN_dummy", subset(missionvote, party == "0")),
          count_votes("EU_dummy", subset(missionvote, party == "0")),
          count_votes("Other_dummy", subset(missionvote, party == "0"))),
  Greens = c(count_votes("NATO_dummy", subset(missionvote, party == "3")),
             count_votes("UN_dummy", subset(missionvote, party == "3")),
             count_votes("EU_dummy", subset(missionvote, party == "3")),
             count_votes("Other_dummy", subset(missionvote, party == "3"))),
  CDU_CSU = c(count_votes("NATO_dummy", subset(missionvote, party == "4")),
              count_votes("UN_dummy", subset(missionvote, party == "4")),
              count_votes("EU_dummy", subset(missionvote, party == "4")),
              count_votes("Other_dummy", subset(missionvote, party == "4"))),
  Left = c(count_votes("NATO_dummy", subset(missionvote, party == "6")),
           count_votes("UN_dummy", subset(missionvote, party == "6")),
           count_votes("EU_dummy", subset(missionvote, party == "6")),
           count_votes("Other_dummy", subset(missionvote, party == "6"))),
  FDP = c(count_votes("NATO_dummy", subset(missionvote, party == "13")),
          count_votes("UN_dummy", subset(missionvote, party == "13")),
          count_votes("EU_dummy", subset(missionvote, party == "13")),
          count_votes("Other_dummy", subset(missionvote, party == "13"))),
  SPD = c(count_votes("NATO_dummy", subset(missionvote, party == "23")),
          count_votes("UN_dummy", subset(missionvote, party == "23")),
          count_votes("EU_dummy", subset(missionvote, party == "23")),
          count_votes("Other_dummy", subset(missionvote, party == "23")))
)

library(knitr)
library(kableExtra)
kable(table_votes, format = "html", caption = "Cumulative Voting Results by Parties and International Organizations") %>%
  kable_styling(full_width = F)


```

```{r}
#Table 1 Main text - Linear model vs. curvilinear model
# Model 1: Linear model with robust standard errors
model1 <- glm(support_dummy ~ bu01o, data = missionvote, family = binomial(link = "logit"))
summary(model1)
coeftest(model1, vcov = vcovHC(model1, type = "HC0"))

coeftest(model1, vcov = vcovCL(model1, cluster = missionvote$party))

# Model 2: Linear model with robust standard errors and control variables
model2 <- glm(support_dummy ~ bu01o + NATO + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile, 
              data = missionvote, family = binomial(link = "logit"))
summary(model2)
coeftest(model2, vcov = vcovHC(model2, type = "HC0"))

coeftest(model2, vcov = vcovCL(model2, cluster = missionvote$party))

# Model 3: Linear model with IO's
model3 <- glm(support_dummy ~ NATO, 
              data = missionvote, family = binomial(link = "logit"))
summary(model3)
coeftest(model3, vcov = vcovHC(model3, type = "HC0"))

coeftest(model3, vcov = vcovCL(model3, cluster = missionvote$party))


# Model 4: Polynomial model with robust standard errors
model4 <- glm(support_dummy ~ poly(bu01o, 2), data = missionvote, family = binomial(link = "logit"))
summary(model4)
# Obtain robust standard errors
coeftest(model4, vcov = vcovHC(model4, type = "HC3"))
coeftest(model4, vcov = vcovCL(model4, cluster = missionvote$party))

# Modell 5: Polynomial model with robust standard errors and Control Variables 
model5 <- glm(support_dummy ~ poly(bu01o, 2) + NATO + avg_govt + + avg_value + avg_unemployment + avg_range_parl_rile, 
              data = missionvote, 
              family = binomial(link = "logit"))
summary(model5)
coeftest(model5, vcov = vcovHC(model5, type = "HC3"))
coeftest(model5, vcov = vcovCL(model5, cluster = missionvote$party))
```

```{r}
#Robustness test CHES
# Model 1: Linear model with robust standard errors
model1 <- glm(support_dummy ~ lrgen, data = missionvote, family = binomial(link = "logit"))
summary(model1)
coeftest(model1, vcov = vcovHC(model1, type = "HC0"))

coeftest(model1, vcov = vcovCL(model1, cluster = missionvote$party))

# Model 2: Linear model with robust standard errors and control variables
model2 <- glm(support_dummy ~ lrgen + NATO + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile, 
              data = missionvote, family = binomial(link = "logit"))
summary(model2)
coeftest(model2, vcov = vcovHC(model2, type = "HC0"))

coeftest(model2, vcov = vcovCL(model2, cluster = missionvote$party))

# Model 3: Linear model with IOs
model3 <- glm(support_dummy ~ NATO, 
              data = missionvote, family = binomial(link = "logit"))
summary(model3)
coeftest(model3, vcov = vcovHC(model3, type = "HC0"))

coeftest(model3, vcov = vcovCL(model3, cluster = missionvote$party))


# Model 4: Polynomial model with robust standard errors
model4 <- glm(support_dummy ~ poly(lrgen, 2), data = missionvote, family = binomial(link = "logit"))
summary(model4)
# Obtain robust standard errors
coeftest(model4, vcov = vcovHC(model4, type = "HC3"))
coeftest(model4, vcov = vcovCL(model4, cluster = missionvote$party))

# Model 5: Polynomial model with robust standard errors and control variables
model5 <- glm(support_dummy ~ poly(lrgen, 2) + NATO + avg_govt + + avg_value + avg_unemployment + avg_range_parl_rile, 
              data = missionvote, 
              family = binomial(link = "logit"))
summary(model5)
coeftest(model5, vcov = vcovHC(model5, type = "HC3"))
coeftest(model5, vcov = vcovCL(model5, cluster = missionvote$party))
```

```{r}
# Graphical representation of Model 4
# Model 4: Polynomial model with robust standard errors n
model4 <- glm(support_dummy ~ poly(bu01o, 2), data = missionvote, family = binomial(link = "logit"))
summary(model4)
# Obtain robust standard errors
coeftest(model4, vcov = vcovHC(model4, type = "HC3"))
coeftest(model4, vcov = vcovCL(model4, cluster = missionvote$party))


new_data <- data.frame(bu01o = seq(min(missionvote$bu01o), max(missionvote$bu01o), length.out = 1000))

predictions <- predict(model4, newdata = new_data, type = "link", se.fit = TRUE)
crit_val <- qnorm(0.975)  
new_data$lower <- predictions$fit - (crit_val * predictions$se.fit)
new_data$upper <- predictions$fit + (crit_val * predictions$se.fit)
new_data$predicted_prob <- plogis(predictions$fit)  
new_data$lower_prob <- plogis(new_data$lower)
new_data$upper_prob <- plogis(new_data$upper)

# Plotting predictions and confidence intervals
p <- ggplot(new_data, aes(x = bu01o, y = predicted_prob)) +
  geom_line(color = "blue", size = 1) +  
  geom_ribbon(aes(ymin = lower_prob, ymax = upper_prob), fill = "blue", alpha = 0.15) +  
  labs(title = "",
       x = "RILE Index",
       y = "Probability (yes-vote = 1)") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0))  

p


ggsave("predicted_probability_plot.png", plot = p, width = 8, height = 6, dpi = 300)

# Robustness test without the AfD table in the appendix + graph
# Data cleaned to exclude AfD
missionvote_no_afd <- subset(missionvote, party != 0)

# Logistic regression model without AfD
model_no_afd <- glm(support_dummy ~ poly(bu01o, 2), data = missionvote_no_afd, family = binomial(link = "logit"))

newdata_no_afd <- data.frame(bu01o = seq(min(missionvote_no_afd$bu01o, na.rm = TRUE), 
                                         max(missionvote_no_afd$bu01o, na.rm = TRUE), 
                                         length.out = 1000))

predictions_no_afd <- predict(model_no_afd, newdata = newdata_no_afd, type = "link", se.fit = TRUE)
crit_val <- qnorm(0.975)  
newdata_no_afd$lower <- predictions_no_afd$fit - (crit_val * predictions_no_afd$se.fit)
newdata_no_afd$upper <- predictions_no_afd$fit + (crit_val * predictions_no_afd$se.fit)
newdata_no_afd$predicted_prob <- plogis(predictions_no_afd$fit)
newdata_no_afd$lower_prob <- plogis(newdata_no_afd$lower)
newdata_no_afd$upper_prob <- plogis(newdata_no_afd$upper)

# Plotting predictions and confidence intervals
ggplot(newdata_no_afd, aes(x = bu01o, y = predicted_prob)) +
  geom_line(color = "blue", size = 1) +  
  geom_ribbon(aes(ymin = lower_prob, ymax = upper_prob), fill = "blue", alpha = 0.15) +  
  labs(title = "",
       x = "RILE Index (Without AfD)",
       y = "Probability (yes-vote = 1)") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0))  

ggsave("C:/Users/Stremming/Documents/Paper 1/Stuff/predicted_probability_no_afd.png", plot = last_plot(), width = 8, height = 6, dpi = 300)

# Table Appendix - without Afd, clustered standard errors and control variables
model_no_afd2 <- glm(support_dummy ~ poly(bu01o, 2) + NATO + avg_govt + + avg_value + avg_unemployment + avg_range_parl_rile, data = missionvote_no_afd, family = binomial(link = "logit"))
summary(model_no_afd2)

coeftest(model_no_afd2, vcov = vcovHC(model_no_afd2, type = "HC3"))
coeftest(model_no_afd2, vcov = vcovCL(model_no_afd2, cluster = missionvote_no_afd$party))

```

```{r} 
#Robustness Tests
#Table in the appendix: Party families with robust standard errors
modelC1 <- glm(support_dummy ~ party_family_PDVD, data = missionvote, family = binomial(link = "logit"))
summary(modelC1)

#Table in the appendix: Robust standard errors for party families and control variables
coeftest(modelC1, vcov = vcovHC(modelC1, type = "HC0"))

modelC2 <- glm(support_dummy ~ party_family_PDVD + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile, 
               data = missionvote, family = binomial(link = "logit"))
summary(modelC2)
coeftest(modelC2, vcov = vcovHC(modelC2, type = "HC0"))

#Table in the appendix: Individual models for IOs
# Conversion of the NATO variable into dummy variables
missionvote$NATO_dummy <- ifelse(missionvote$NATO == 1, 1, 0)  # NATO Dummy
missionvote$UN_dummy <- ifelse(missionvote$NATO == 2, 1, 0)    # UN Dummy
missionvote$EU_dummy <- ifelse(missionvote$NATO == 3, 1, 0)    # EU Dummy
missionvote$Other_dummy <- ifelse(missionvote$NATO == 8, 1, 0) # Others Dummy

#NATO Mission
model_nato <- glm(support_dummy ~ bu01o + avg_value +  avg_govt + avg_unemployment + avg_range_parl_rile,
                data = subset(missionvote, UN_dummy == 1), 
                family = binomial(link = "logit"))
summary(model_nato)
clustered_se_nato <- vcovCL(model_nato, cluster = ~ party)
coeftest(model_nato, vcov = clustered_se_nato)
model_nato <- glm(support_dummy ~ poly(bu01o, 2) + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile,
                  data = subset(missionvote, NATO_dummy == 1), 
                  family = binomial(link = "logit"))
summary(model_nato)

clustered_se_nato <- vcovCL(model_nato, cluster = ~ party)
coeftest(model_nato, vcov = clustered_se_nato)

# UN Mission
model_un <- glm(support_dummy ~ poly(bu01o, 2) +avg_value +  avg_govt + avg_unemployment + avg_range_parl_rile,
                data = subset(missionvote, UN_dummy == 1), 
                family = binomial(link = "logit"))
summary(model_un)
clustered_se_un <- vcovCL(model_un, cluster = ~ party)
coeftest(model_un, vcov = clustered_se_un)

# EU Mission
model_eu <- glm(support_dummy ~ poly(bu01o, 2) + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile,
                data = subset(missionvote, EU_dummy == 1), 
                family = binomial(link = "logit"))
summary(model_eu)
clustered_se_eu <- vcovCL(model_eu, cluster = ~ party)
coeftest(model_eu, vcov = clustered_se_eu)


# Table in the appendix: Interaction between bu01o and the dummy variables for NATO, UN and EU
model_interaction1 <- glm(support_dummy ~ bu01o * NATO_dummy + bu01o * UN_dummy + bu01o * EU_dummy, 
                          data = missionvote, family = binomial(link = "logit"))
summary(model_interaction1)
clustered_se <- vcovCL(model_interaction1, cluster = ~ party)
coeftest(model_interaction1, vcov = clustered_se)


model_interaction2 <- glm(support_dummy ~ bu01o * NATO_dummy + bu01o * UN_dummy + bu01o * EU_dummy +  avg_value + avg_govt + avg_unemployment + avg_range_parl_rile, 
                         data = missionvote, family = binomial(link = "logit"))
summary(model_interaction2)

clustered_se <- vcovCL(model_interaction2, cluster = ~ party)

coeftest(model_interaction2, vcov = clustered_se)

model_interaction_poly <- glm(support_dummy ~ poly(bu01o, 2) * NATO_dummy + 
                              poly(bu01o, 2) * UN_dummy + 
                              poly(bu01o, 2) * EU_dummy + 
                              avg_govt + avg_unemployment + avg_range_parl_rile, 
                              data = missionvote, family = binomial(link = "logit"))
summary(model_interaction_poly)
clustered_se_poly <- vcovCL(model_interaction_poly, cluster = ~ party)
coeftest(model_interaction_poly, vcov = clustered_se_poly)

model_nato <- glm(support_dummy ~ poly(bu01o, 2) * NATO_dummy + avg_value + 
                  avg_govt + avg_unemployment + avg_range_parl_rile, 
                  data = missionvote, family = binomial(link = "logit"))
clustered_se_nato <- vcovCL(model_nato, cluster = ~ party)
coeftest(model_nato, vcov = clustered_se_nato)

model_un <- glm(support_dummy ~ poly(bu01o, 2) * UN_dummy + avg_value + 
                avg_govt + avg_unemployment + avg_range_parl_rile, 
                data = missionvote, family = binomial(link = "logit"))
summary(model_un)
clustered_se_un <- vcovCL(model_un, cluster = ~ party)
coeftest(model_un, vcov = clustered_se_un)

model_eu <- glm(support_dummy ~ poly(bu01o, 2) * EU_dummy + avg_value + 
                avg_govt + avg_unemployment + avg_range_parl_rile, 
                data = missionvote, family = binomial(link = "logit"))
summary(model_eu)
clustered_se_eu <- vcovCL(model_eu, cluster = ~ party)
coeftest(model_eu, vcov = clustered_se_eu)

```

```{r} 
#Additional Visualisation
missionvote$org <- NA
missionvote$org[missionvote$NATO_dummy == 1] <- "NATO"
missionvote$org[missionvote$UN_dummy == 1] <- "UN"
missionvote$org[missionvote$EU_dummy == 1] <- "EU"
missionvote$org[missionvote$Other_dummy == 1] <- "Others"

ggplot(missionvote, aes(x = bu01o, y = support_dummy, color = org)) +
  geom_point(alpha = 0.7) +
  facet_wrap(~ org) +
  labs(title = "",
       x = "Ideological Position (RILE Index)",
       y = "Probability (Mission Vote)",
       color = "International Organizations") +
  scale_color_manual(values = c("blue", "green", "red", "purple")) +  
  theme_minimal()

ggsave("C:/Users/Stremming/Documents/Paper 1/Stuff/ideological_position_mission_vote.png", width = 10, height = 7, dpi = 300)

model_interaction_cat <- glm(support_dummy ~ poly(bu01o, 2) * NATO, data = missionvote, family = binomial(link = "logit"))

missionvote$predicted_prob_cat <- predict(model_interaction_cat, type = "response")

plot_data <- missionvote
plot_data$NATO <- as.factor(plot_data$NATO) 
ggplot(plot_data, aes(x = bu01o, y = predicted_prob_cat, color = NATO)) +
  geom_line(size = 1) +  
  labs(
    title = "Interaction between Ideology and NATO on support",
    x = "Ideology",
    y = "Predicted probability"
  ) +
  theme_minimal()


new_data <- missionvote %>%
  group_by(NATO) %>%
  do(data.frame(bu01o = seq(min(missionvote$bu01o), max(missionvote$bu01o), length.out = 100))) %>%
  ungroup()

conf_int <- predict(model_interaction_cat, newdata = new_data, type = "link", se.fit = TRUE)

new_data <- new_data %>%
  mutate(predicted_prob = plogis(conf_int$fit), 
         lower = plogis(conf_int$fit - 1.96 * conf_int$se.fit),  
         upper = plogis(conf_int$fit + 1.96 * conf_int$se.fit))  

ggplot(new_data, aes(x = bu01o, y = predicted_prob, color = factor(NATO))) +
  geom_line(size = 1) +  
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = factor(NATO)), alpha = 0.1, color = NA) + 
  labs(
    title = "",
    x = "Ideological Position",
    y = "Probability Mission Vote",
    color = "International Organizations",  
    fill = "International Organizations"     #
  ) +
  scale_color_manual(
    values = c("red", "green", "blue", "purple"),
    labels = c("8" = "Others", "1" = "NATO", "2" = "UN", "3" = "EU")  
  ) +
  scale_fill_manual(
    values = c("red", "green", "blue", "purple"),
    labels = c("8" = "Others", "1" = "NATO", "2" = "UN", "3" = "EU")  
  ) +
  theme_minimal()

ggsave("C:/Users/Stremming/Documents/Paper 1/Stuff/interaction.eps", width = 10, height = 7, dpi = 300)
```


```{r}
# Additional Models
model_interaction_voll <- glm(support_dummy ~ poly(bu01o, 2) * NATO + avg_govt + avg_value + avg_unemployment + avg_range_parl_rile, 
                              data = missionvote, family = binomial(link = "logit"))

# Calculate coefficients and confidence intervals
coef_data <- tidy(model_interaction_voll, conf.int = TRUE, conf.level = 0.95)

# Calculate robust standard errors
coef_data$std.error <- sqrt(diag(vcovHC(model_interaction_voll, type = "HC3")))
coef_data$conf.low <- coef_data$estimate - 1.96 * coef_data$std.error
coef_data$conf.high <- coef_data$estimate + 1.96 * coef_data$std.error

term_labels <- c(
  "poly(bu01o, 2)2:NATO3" = "RILE^2*EU",
  "poly(bu01o, 2)1:NATO3" = "RILE*EU",
  "poly(bu01o, 2)2:NATO2" = "RILE^2*UN",
  "poly(bu01o, 2)1:NATO2" = "RILE*UN",
  "poly(bu01o, 2)2:NATO1" = "RILE^2*NATO",
  "poly(bu01o, 2)1:NATO1" = "RILE*NATO",
  "avg_range_parl_rile" = "Range of Parties",
  "avg_unemployment" = "Unemployment",
  "avg_value" = "Party Position",
  "avg_govt" = "Government",
  "NATO3" = "EU",
  "NATO2" = "UN",
  "NATO1" = "NATO",
  "poly(bu01o, 2)2" = "RILE^2",
  "poly(bu01o, 2)1" = "RILE",  # Linearer Term für RILE
  "(Intercept)" = "Intercept"
)

coef_data$term <- term_labels[coef_data$term]

print("Terms after renaming:")
print(unique(coef_data$term))
coef_data$term[is.na(coef_data$term)] <- "Unlabeled"
term_order <- c(
  "Intercept", "Range of Parties", "Unemployment", "Party Position", 
  "Government", "RILE*EU", "RILE^2*EU", "RILE^2*UN", "RILE^2*NATO", 
  "RILE*UN", "RILE*NATO", "EU", "UN", "NATO", "RILE^2", "RILE"
)

coef_data$term <- factor(coef_data$term, levels = term_order)

ggplot(coef_data, aes(x = estimate, y = term)) +
  geom_point(size = 3, color = "darkgray") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2, color = "darkgray") +
  labs(
    title = "",
    x = "Estimate",
    y = "Terms"
  ) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 10))  # Schriftgröße der Achsentitel anpassen

ggsave("C:/Users/Stremming/Documents/Paper 1/Stuff/interaction_plot.png", width = 10, height = 7, dpi = 300)

```

```{r}
# Model for the period 1990-2013
model5_1990_2013 <- glm(support_dummy ~ poly(bu01o, 2) + NATO + avg_govt + avg_value + avg_unemployment + avg_range_parl_rile, 
                        data = subset(missionvote, date >= 1990 & date <= 2013), 
                        family = binomial(link = "logit"))

# Clustered standard errors by party for 1990-2013
cluster_1990_2013 <- missionvote$party[missionvote$date >= 1990 & missionvote$date <= 2013]
coeftest(model5_1990_2013, vcov = vcovCL(model5_1990_2013, cluster = cluster_1990_2013))


# Model for the period 2014-2019
model5_2014_2019 <- glm(support_dummy ~ poly(bu01o, 2) + NATO + avg_govt + avg_value + avg_unemployment + avg_range_parl_rile, 
                        data = subset(missionvote, date >= 2014 & date <= 2019), 
                        family = binomial(link = "logit"))

# Clustered standard errors by party for 2014-2019
cluster_2014_2019 <- missionvote$party[missionvote$date >= 2014 & missionvote$date <= 2019]
coeftest(model5_2014_2019, vcov = vcovCL(model5_2014_2019, cluster = cluster_2014_2019))

```

```{r}

# Calculate Cook's Distance
cooks_distance <- cooks.distance(model5)
plot(cooks_distance, pch = 20, main = "Cook's Distance Plot", ylab = "Cook's Distance")
abline(h = 4 / nrow(missionvote), col = "red")  
text(x = 1:length(cooks_distance), y = cooks_distance, labels = ifelse(cooks_distance > 4/nrow(missionvote), names(cooks_distance), ""), pos = 4, cex = 0.7, col = "red")

# Calculate leverage and standardized residuals
leverage <- hatvalues(model5)
standardized_residuals <- rstandard(model_interaction_voll)

## Scatterplot of leverage versus standardized residuals
plot(leverage, standardized_residuals, 
     main = "Leverage vs. Standardized Residuals", 
     xlab = "Leverage", ylab = "Standardized Residuals", pch = 20)
abline(h = 0, col = "red", lty = 2)  
abline(v = 2 * mean(leverage), col = "blue", lty = 2)  

# Identification of influential points
text(x = leverage, y = standardized_residuals, labels = ifelse(leverage > 2 * mean(leverage), names(leverage), ""), pos = 4, cex = 0.7, col = "blue")


# Calculate Cook's Distance
cooks_distance <- cooks.distance(model_interaction_voll)
plot(cooks_distance, pch = 50, main = "Cook's Distance Plot", ylab = "Cook's Distance")
abline(h = 4 / nrow(missionvote), col = "red")  
text(x = 1:length(cooks_distance), y = cooks_distance, labels = ifelse(cooks_distance > 4/nrow(missionvote), names(cooks_distance), ""), pos = 4, cex = 0.7, col = "red")

# Calculate leverage and standardized residuals
leverage <- hatvalues(model_interaction_voll)
standardized_residuals <- rstandard(model_interaction_voll)

# Scatterplot of leverage versus standardized residuals
plot(leverage, standardized_residuals, 
     main = "Leverage vs. Standardized Residuals", 
     xlab = "Leverage", ylab = "Standardized Residuals", pch = 20)
abline(h = 0, col = "red", lty = 2)  
abline(v = 2 * mean(leverage), col = "blue", lty = 2)  

# Identification of influential points
text(x = leverage, y = standardized_residuals, labels = ifelse(leverage > 2 * mean(leverage), names(leverage), ""), pos = 4, cex = 0.7, col = "blue")

```

```{r}
#DV alternative modeling
#Model 1: Linear model with robust standard errors
model1 <- glm(support_dummy1 ~ bu01o, data = missionvote, family = binomial(link = "logit"))
summary(model1)
coeftest(model1, vcov = vcovHC(model1, type = "HC0"))

coeftest(model1, vcov = vcovCL(model1, cluster = missionvote$party))

#Model 2: Linear model with robust standard errors and control variables
model2 <- glm(support_dummy1 ~ bu01o + NATO + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile, 
              data = missionvote, family = binomial(link = "logit"))
summary(model2)
coeftest(model2, vcov = vcovHC(model2, type = "HC0"))

coeftest(model2, vcov = vcovCL(model2, cluster = missionvote$party))

# Model 3: Linear model with IOs
model3 <- glm(support_dummy1 ~ NATO  + avg_value + avg_govt + avg_unemployment + avg_range_parl_rile,
              data = missionvote, family = binomial(link = "logit"))
summary(model3)
coeftest(model3, vcov = vcovHC(model3, type = "HC0"))

coeftest(model3, vcov = vcovCL(model3, cluster = missionvote$party))


# Model 4: Polynomial model with robust standard errors
model4 <- glm(support_dummy1 ~ poly(bu01o, 2), data = missionvote, family = binomial(link = "logit"))
summary(model4)
# Obtain robust standard errors
coeftest(model4, vcov = vcovHC(model4, type = "HC3"))
coeftest(model4, vcov = vcovCL(model4, cluster = missionvote$party))

# Model 5: Polynomial model with robust standard errors and control variables
model5 <- glm(support_dummy1 ~ poly(bu01o, 2) + NATO + avg_govt + + avg_value + avg_unemployment + avg_range_parl_rile, 
              data = missionvote, 
              family = binomial(link = "logit"))
summary(model5)
coeftest(model5, vcov = vcovHC(model5, type = "HC3"))
coeftest(model5, vcov = vcovCL(model5, cluster = missionvote$party))


```